function epsc_mash(varargin)
% Calculate the amplitude of the EPSC based on the prior EPSCs
% 

sf = getmainselection; % call the routine that lists the selected record indices
if(sf > 0) % make sure there's a selection
    QueMessage('EPSC_mash', 1); % clear the que and show a message
    for i = 1:length(sf) % loop through the selected database entries
        epsc_mash2(sf(i)); % call the analysis routine with argumetns
    end;
end;

return;

function epsc_mash2(sf)
global ALLCH

[s d] = block_info2(sf);
SR=d.Sample_Rate.v*2;
nrec = size(ALLCH{2}, 1);
d.rate=(SR)*ones(nrec, 1);
d.nr_channel = d.Channels.v(1);
d.mode = 5;
d.record = 1:nrec;
d.frec = 1;
d.lrec = nrec;
d.nr_points = d.Points.v;

% other parameters
tmax = 100; % max recovery window to watch
t_ssstart = 100; % time at which steady-state rate is reached

t=0:0.1:tmax;
pevent = 2;
poststart = 500;
minpost = 0;

[fs, fisi, nrs, sptrain] = find_spikes2(d, ALLCH{2}, 0, poststart, [4000, -1]);
[fsr, fisir, nrsr, sptrainr] = find_spikes2(d, ALLCH{2}, minpost+poststart, 1050, [4000, -1]);

for i = 1:length(sptrainr)
    sptrainr(i).latency
end;



[hpanel] = multiplot('epsc_recovery', 2, 3, 1, [0.05 0.05 0.9 0.9], [5 5]);

%newfigure('raster', 'Raster');
axes(hpanel(2));
rv = cell(nrec, 1);
ev = cell(nrec, 1);
for i = 1:nrec
    sl = [sptrain(i).latency];
    sl = sl(find(~isnan(sl)));
    rv{i} = i+zeros(length(sl), 1);
    ev{i} = sl;
    plot([ev{i}], [rv{i}], 'r.');
    hold on;
end;

%%
tw = floor(0.8/(SR/1000)); % total window
dw = floor(0.5/(SR/1000)); % deadwindow


amp = cell(nrec, 1);
wv=amp;
for i = 1:nrec
    bl = mean(ALLCH{2}(i, 1:25));
    for j = 1:length(ev{i})
        it1 = floor(ev{i}(j)/(SR/1000))+dw;
        it2 = it1 + tw;
        amp{i}(j) = min(ALLCH{2}(i, it1:it2)-bl);
        wv{i,j} = ALLCH{2}(i, it1:it2)-bl;
    end;
end;
%
% now calculate normalized amplitudes. Use the mean of the first response
% across all trials for the normalization
%
amp_first = zeros(nrec, 1);
for i = 1:nrec;
    amp_first(i) = amp{i}(1);
end
amp_mean = mean(amp_first);

% this calculates the post-stimulus recovery tests, if they are present
amp_post = cell(nrec, 1);  % amplitudes
wv_post = amp; % waveforms
lat_post = amp;
for i = 1:nrec
    bl = mean(ALLCH{2}(i, 1:25));
    if(sptrainr(i).source == i)
        lat = sptrainr(i).latency;
        lat = lat(1)+minpost+poststart; % only the first latency.
        it1 = floor(lat/(SR/1000))+dw;
        it2 = it1 + tw;
        amp_post{i} = min(ALLCH{2}(i, it1:it2)-bl)/amp_mean;
        lat_post{i} = lat-poststart;
        wv{i} = ALLCH{2}(i, it1:it2)-bl;
    end;
end;



%newfigure('amp', 'Amplitudes');
axes(hpanel(1));
for i = 1:nrec
    for j = 1:length(wv{i})
        plot(wv{i,j});
        hold on;
    end;
end;

axes(hpanel(3));
for j = 1:nrec
    plot([ev{j}], [amp{j}]/amp_mean, 'k.');
    hold on;
end;

u = get(gca, 'Ylim');
set(gca, 'Ylim', [0 max(u)]);
%%
% events versus prior isi...
%
axes(hpanel(4));
%newfigure('priors', 'Prior ISI');
sym={'r.', 'y.', 'g.', 'c.', 'b.', 'm.', 'rx', 'yx', 'gx', 'cx', 'bx', 'mx', 'rs', 'ys', 'gs', 'cs', 'bx', 'ms'};

lin={'r-', 'y-', 'g-', 'c-', 'b-', 'm-', 'r--', 'y--', 'g--', 'c--', 'b--', 'm--', 'r:', 'y:', 'g:', 'c:', 'b:', 'm:'};

allev=cell(max(pevent), 1);
allet= allev;
for i = 1:nrec
    evt = [ev{i}]; % get event times for this run
    it=find(evt > t_ssstart); % eliminate non-ss part of response
    evt = evt(it);
    eva = amp{i}/amp_mean; % get amplitudes for this run
    eva = eva(it);
    evtd = diff(evt); % prior ISIs
    plot(evtd, eva(2:end), '.k'); % amplitdues are in yellow
    hold on
    evtd2=[];
    eva2=[];
    for k = pevent % prior event count...
        for j = k:length(eva)
            evtd2(j-k+1) = mean(diff(evt(j-k+1:j))); %#ok<AGROW>
            eva2(j-k+1) = eva(j); %#ok<AGROW>
        end;
        allet{k} = [allet{k} evtd2];
        allev{k} = [allev{k} eva2];
    end;

end;

maxiter = 10000;

for k = pevent
    X = sortrows([allet{k}; allev{k}]');
    a0 = 1/3; a1 = -2/3; tau = 25;
    [fpar] = mrqfit('exponential', [a0 a1 tau], X(:,1), X(:,2), [], ...
        [], [-2 -2 1], [2 2 2*max(t)], maxiter, []); % VP LP UB Imax Tol

        a0 = fpar(1); a1=fpar(2); tau = fpar(3);
        plot(X(:,1), X(:,2), sym{k});
        yfit = a0 + a1*exp(-t/tau);
        plot(t, yfit, lin{k});
 
end;


%newfigure('ISIGroupdT', 'Recovery grouped by prior ISI(n)');
axes(hpanel(5));
mean_isi = zeros(nrec, 1);
priordt = mean_isi;
thisamp = mean_isi;
k = 1;
for i  = 1:nrec
    evt = [ev{i}];
    for j = 1:length(evt);
        if(evt(j) < t_ssstart)
            continue;
        end;
        if(j <= pevent+1)
            continue;
        end;
        mean_isi(k) = mean(diff(evt(j-pevent-1:j-1))); % mean isi of prior pevent events
        priordt(k) = evt(j)-evt(j-1);
        thisamp(k) = amp{i}(j)/amp_mean;
        k = k + 1;
    end;
end;
vlist = [4 8 12 16 20];
tfit = 0:0.1:50;
nlist = length(vlist)-1;
tau = 25*ones(nlist, 1);
a1 = 1/3*ones(nlist, 1);
a0 = 2/3*ones(nlist, 1);
v = cell(nlist, 1);
for i = 1:nlist
    v{i} = find(mean_isi > vlist(i) & mean_isi <= vlist(i+1));
    if(isempty(v{i}) || length(v{i}) < 3 || (max(v{i})-min(v{i})) < 15)
        continue;
    end;
    
        plot(priordt(v{i}), thisamp(v{i}), sym{i});
    hold on
    % check out a fit
    a0(i) = 1/3; a1(i) = -2/3; tau(i) = 25;
    [fpar, chisq, niter] = mrqfit('exponential', [a0(i) a1(i) tau(i)], priordt(v{i}), thisamp(v{i}), [], ...
        [], [-2 -2 2*min(priordt(v{i}))], [2 2 2*max(t)], 10000, []); % VP LP UB Imax Tol
    a0(i) = fpar(1);
    a1(i) = fpar(2);
    tau(i) = fpar(3);
    %plot((:,1), X(:,2), sym{k});
    yfit = a0(i) + a1(i)*exp(-tfit/tau(i));
    plot(tfit, yfit, lin{i});
    for j = 1:length(fpar)
        fprintf(1, 'fpar(%d): %8.3f \n', i, fpar(j));
    end;
    fprintf(1, 'chisq: %f  niter = %d\n\n', chisq, niter);
    [m, b, r, p] = linreg(priordt(v{i}), thisamp(v{i}));
    mse_exp = 0; mse_lin = 0;
    for j = 1:length(thisamp(v{i}))
        mse_exp = mse_exp + (thisamp(v{i}(j)) - (a0(i)+a1(i)*exp(-priordt(v{i}(j))/tau(i))))^2;
        mse_lin = mse_lin + (thisamp(v{i}(j)) - (m * priordt(v{i}(j)) + b))^2;
    end;
    fprintf(1, 'u(int) = %f  Exp Err: %f    lin Err: %f\n', mean(vlist(i:i+1)), mse_exp, mse_lin);
end;
%
% show the time constant of the fit...
axes(hpanel(6));
latp = cell2mat(lat_post);
latd = diff(latp);
ampp = cell2mat(amp_post);
isk = find(latd > 1);
isk = [isk; length(latp)];
k = 1;
delay = zeros(length(isk), 1);
irec = delay;
sirec = delay;
for i = 1:length(isk)
    delay(i) = mean(latp(k:isk(i)));
    [irec(i) sirec(i)] = mean_var(ampp(k:isk(i)));
    k = k + isk(1);
end;

plot_witherror(delay, irec, sqrt(sirec));

%for i = 1:nlist
%    if(~isempty(v{i}))
%        plot(mean(mean_isi(v{i})), tau(i), 'ko');
%    end;
%    hold on
%end;



